Ejercicios del libro de Faraway

1. (Ejercicio 1 cap. 6 pág. 97)

Using the sat dataset, fit a model with the total SAT score as the response and expend, salary, ratio and takers as predictors. Perform regression diagnostics on this model to answer the following questions. Display any plots that are relevant. Do not provide any plots about which you have nothing to say. Suggest possible improvements or corrections to the model where appropriate.

n <- dim(sat)[1]
p <- dim(sat)[2]
model <- 
    lm(total ~ expend + salary + ratio + takers, data = sat)

sat %>%
    select(total, expend, salary, ratio, takers) %>%
    ggpairs(progress = F)

Lo primero que hago es crear una matriz de dispersión para tener una idea de como se distribuyen las variables y la relación 2x2 que tienen entre si. Utilizo la librería GGally que tiene la función ggpairs.

  1. Check the constant variance assumption for the errors.
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos

Gráficos Diagnósticos

Para evaluar la homocedasticidad del error utilizo los gráficos de:

  • residuales v/s valor ajustado
  • residuales estandarisados v/s valor ajustado.

En este caso se observa cierta irregularidad en el gráfico \(\hat{y} \times \hat{\epsilon}\) que podría traducir no-linearidad, aunque dentro de todo no se observa un patrón definido, ni tan anormal que haga pensar que el erro no sigue una homocedasticidad.

  1. Check the normality assumption.
shapiro.test(resid(model)) %>%
  pander()
Shapiro-Wilk normality test: resid(model)
Test statistic P value
0.9769 0.4304

Para evaluar la normalidad del error, utilizo el gráfico de Q-Q de los residuales. En este caso el modelo se acerca bastante a la linea recta, con algunos valores raros que se alejan, que corresponden a los con valores mas altos de residuales, y que habría que examinar mejor.

Por otra parte podemos utilizar el test de Shapiro-Wilk cuya H0 es que el error no difiere de la distribución normal. En este caso no hay evidencia para rechazar la H0

  1. Check for large leverage points.
halfnorm(hatvalues(model), labs = rownames(sat))

augment(model) %>%
    filter(.hat > 2*(p/n)) %>%
  pander()
.rownames total expend salary ratio takers .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
California 902 4.992 41.08 24 45 917.8 17.37 -15.85 0.2821 32.95 0.02571 -0.572
Utah 1076 3.656 29.08 24.3 4 1010 17.67 65.77 0.2921 30.9 0.4715 2.39

Para evaluar los valores extremos en el modelo, utilizamos dos gráficos.

  • La gráfica de los valores ajustados versus los residuales estandarizados (\(\hat{y} \times \sqrt{|{\hat{\epsilon}}|}\))
  • Una gráfica de media normal (halfnormal) versus los residuales estandarizados. De esta manera se pueden identificar los valores mas extremos en \(X\).

En este caso podemos ver que los valores para California y Utah cumplen con las características para valores extremos.

  1. Check for outliers.
(outliers <- 
    augment(model) %>%
    mutate(.t.resid = rstudent(model)) %>%
    arrange(-abs(.t.resid)) %>%
    head(5)) %>%
    pander()
Table continues below
.rownames total expend salary ratio takers .fitted .se.fit .resid .hat .sigma .cooksd
West Virginia 932 6.107 31.94 14.8 17 1023 8.147 -90.53 0.06207 29.92 0.1081
Utah 1076 3.656 29.08 24.3 4 1010 17.67 65.77 0.2921 30.9 0.4715
North Dakota 1107 4.775 26.33 15.3 5 1040 9.31 66.57 0.08105 31.37 0.07954
New Hampshire 935 5.859 34.72 15.6 70 869.1 9.412 65.87 0.08284 31.4 0.07989
Nevada 917 5.16 34.84 18.7 30 971.1 6.967 -54.15 0.04539 32 0.02731
.std.resid .t.resid
-2.859 -3.124
2.39 2.53
2.124 2.214
2.103 2.19
-1.695 -1.732 
abs(outliers$.t.resid) > abs(qt(.05/(n * 2), n - p))
## [1] FALSE FALSE FALSE FALSE FALSE

Para valorar posibles outliers podemos utilizar la estrategia planteada en el libro de Faraway, y calcular los residuales studentizados. En este sentido los residuales cuyo valor superen un limite determinado por la corrección de Bonferroni, se podrían considerar posibles outliers.

En este caso el valor que mas puede ser un outlier es West Virginia, aunque el valor de .t.resid de west virgina no es mayor al valor limite de p corregido.

model_1 <- 
    sat %>%
    rownames_to_column() %>%
    filter(rowname != "West Virginia") %>%
    lm(total ~ expend + salary + ratio + takers, data = .)

pander(model)
Fitting linear model: total ~ expend + salary + ratio + takers
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 1046 52.87 19.78 7.858e-24
expend 4.463 10.55 0.4231 0.6742
salary 1.638 2.387 0.6861 0.4962
ratio -3.624 3.215 -1.127 0.2657
takers -2.904 0.2313 -12.56 2.607e-16 
pander(model_1)
Fitting linear model: total ~ expend + salary + ratio + takers
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 1058 48.51 21.8 3.259e-25
expend 7.359 9.693 0.7592 0.4518
salary 1.088 2.191 0.4963 0.6221
ratio -3.94 2.943 -1.338 0.1876
takers -2.972 0.2127 -13.98 8.685e-18

Como West Virgina es un valor extremo (esta lejos de la linea del 0 en el plot de \(\hat{y} \times \sqrt{|{\hat{\epsilon}}|}\))

  1. Check for influential points.
plot(model, 4)

halfnorm(cooks.distance(model), labs = rownames(sat))

model_2 <- 
    sat %>%
    rownames_to_column() %>%
    filter(rowname != "Utah") %>%
    lm(total ~ expend + salary + ratio + takers, data = .)


pander(model)
Fitting linear model: total ~ expend + salary + ratio + takers
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 1046 52.87 19.78 7.858e-24
expend 4.463 10.55 0.4231 0.6742
salary 1.638 2.387 0.6861 0.4962
ratio -3.624 3.215 -1.127 0.2657
takers -2.904 0.2313 -12.56 2.607e-16 
pander(model_1)
Fitting linear model: total ~ expend + salary + ratio + takers
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 1058 48.51 21.8 3.259e-25
expend 7.359 9.693 0.7592 0.4518
salary 1.088 2.191 0.4963 0.6221
ratio -3.94 2.943 -1.338 0.1876
takers -2.972 0.2127 -13.98 8.685e-18 
pander(model_2)
Fitting linear model: total ~ expend + salary + ratio + takers
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 1094 53.42 20.48 4.043e-24
expend -0.9427 10.19 -0.0925 0.9267
salary 3.096 2.328 1.33 0.1904
ratio -7.639 3.428 -2.229 0.031
takers -2.931 0.2188 -13.4 3.946e-17

Para evaluar los puntos mas influyentes podemos utilizar la distancia de Cook que refleja los cambios en el modelo al no incluir una observación determinada. En este caso destacan Utah y West Virginia como observaciones influyentes. Además podemos volver a mirar el gráfico de Residuales v/s leverage, que también refleja puntos con un grado alto de palanca y un residual alto.

influence(model)$coefficients %>%
    as_tibble(rownames = "id") %>%
    gather(key, value, -id) %>%
    ggplot(aes(fct_rev(id), value, group = key)) +
    geom_line(alpha = .75)+
    coord_flip() +
    geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
    facet_grid(~ key, scales = "free") +
    scale_colour_viridis_d(end=.85)

La función influence permite obtener los coeficientes al restar una observación \(x_i\) (“leave one out”). De esta manera puedo construir un gráfico que revela los cambios en los diferentes coeficientes para cada una de las observaciones. En este gráfico se ve claramente como Utah genera los mayores cambios en todos los coeficientes.

plot(model, 5)

pander(sat[44,])
  expend ratio salary takers verbal math total
Utah 3.656 24.3 29.08 4 513 563 1076

El modelo multivariable es explicativo, con un valor de R^2 relativamente alto. En el modelo el gasto por estudiante expend y el sueldo salary se relacionan de manera positiva con el resultado del test total. La ratio de profesores y el porcentaje de estudiantes que hacen el examen de manera negativa.

Utah es un caso raro en la medida que tiene un score muy alta en relación a un gasto expend y un sueldo salary inferior al promedio por una parte, y una ratio. El porcentaje de alumnos elegibles takers llama la atención, ya que parece extremadamente bajo en relación a los otros parámetros.

  1. Check the structure of the relationship between the predictors and the response.

En primer lugar me remito al gráfico de dispersión de (a), que permite reflejar la distribución y la relacion 1 a 1 de las variables.

# added variable 
car::avPlots(model)
## Registered S3 methods overwritten by 'car':
##   method                          from
##   influence.merMod                lme4
##   cooks.distance.influence.merMod lme4
##   dfbeta.influence.merMod         lme4
##   dfbetas.influence.merMod        lme4

Para evaluar la estructura del modelo son de utilidad los gráficos de regresión parcial. La regresión parcial permite aislar la influencia de una variable independiente con la variable dependiente, restando la influencia de las otras variables independientes. Además permiten valorar la presencia de outliers y puntos influyentes.

En este caso podemos ver la clara relación inversa de la variable takers con el total y la relación mas débil de las otras variables. POdemos observar ademas el caso raro de Utah

par(mfrow = c(2,2))
termplot(model, partial.resid = T)

Los gráficos de residuales parciales, permiten conocer la respuesta aislando el efecto del resto de variables independientes. También permiten poder observar la presencia de diferencias estructurales, de no -linearidad

2. (Ejercicio 2 cap. 6 pág. 97)

Using the teengamb dataset, fit a model with gamble as the response and the other variables as predictors. Answer the questions posed in the previous question.

n <- dim(teengamb)[1]
p <- dim(teengamb)[2]

teengamb_cl <- 
    teengamb %>%
    rownames_to_column(".rownames") %>%
    mutate(.rownames = factor(.rownames), 
           sex = factor(sex, levels = c (0,1), labels = c('male', 'female')))

model <- 
    lm(gamble ~ sex + status + income + verbal, data = teengamb_cl)

teengamb_cl %>%
    select(gamble, sex, status, income, verbal) %>%
    ggpairs(progress = F, aes(color = sex), lower = list(combo = wrap("facethist", bins = 10)))

En este caso gráfico con utilizando sex como un factor.

  1. Check the constant variance assumption for the errors.
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos

Gráficos Diagnósticos

En este caso el gráfico podría definir cierta heterocedasticidad en la medida que a mayor valor de \(\hat{y}\) parece haber mayor dispersión del error.

var.test(resid(model)[teengamb_cl$sex == "female"], resid(model)[teengamb_cl$sex == "male"]) %>%
  pander()
F test to compare two variances: resid(model)[teengamb_cl$sex == "female"] and resid(model)[teengamb_cl$sex == "male"]
Test statistic num df denom df P value Alternative hypothesis ratio of variances
0.316 18 27 0.01374 * two.sided 0.316

En este caso al menos en cuanto a las poblaciones según sexo, hay una diferencia significativa en las varianzas

  1. Check the normality assumption.
shapiro.test(resid(model)) %>%
  pander()
Shapiro-Wilk normality test: resid(model)
Test statistic P value
0.8684 8.16e-05 * * *

En este caso podemos observar que el gráfico Q-Q difiere de la normalidad con un patrón de cola larga. Por otra parte la preba de Shapiro Wilk refleja que se rechaza la H0 de normalidad de la distribución del error.

  1. Check for large leverage points.
halfnorm(hatvalues(model), labs = teengamb_cl$.rownames)

augment(model) %>% 
    mutate(.rownames = teengamb_cl$.rownames) %>%
    filter(.hat > 2*(p/n)) %>%
  pander()
gamble sex status income verbal .fitted .se.fit .resid .hat .sigma .cooksd .std.resid .rownames
88 male 18 12 2 77.12 11.1 10.88 0.2395 22.88 0.01904 0.5498 31
90 male 38 15 7 78.25 10.68 11.75 0.2213 22.87 0.01957 0.5867 33
14.1 male 28 1.5 1 28.5 12.67 -14.4 0.3118 22.8 0.05304 -0.7651 35
69.7 male 61 15 9 73.54 12.46 -3.836 0.3016 22.95 0.003535 -0.2023 42

En este caso podemos ver que los valores para 42 y 35 cumplen con las características para valores extremos. Habría que revisar también los 31 y 33.

  1. Check for outliers.
augment(model) %>%
    mutate(.t.resid = rstudent(model)) %>%
    filter(abs(.t.resid) > abs(qt(.05/(n * 2), n - p))) %>%
    pander()
gamble sex status income verbal .fitted .se.fit .resid .hat .sigma .cooksd .std.resid .t.resid
156 male 27 10 4 61.75 7.984 94.25 0.1238 16.74 0.5565 4.438 6.016

En este caso el valor 24 parece ser un outlier

model_1 <- 
    teengamb_cl %>%
    filter(.rownames != "24") %>%
    lm(gamble ~ sex + status + income + verbal, data = .) %>%
    pander()

pander(model)
Fitting linear model: gamble ~ sex + status + income + verbal
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.56 17.2 1.312 0.1968
sexfemale -22.12 8.211 -2.694 0.01011
status 0.05223 0.2811 0.1858 0.8535
income 4.962 1.025 4.839 1.792e-05
verbal -2.959 2.172 -1.362 0.1803 
pander(model_1)
## Warning in pander.default(model_1): No pander.method for "knit_asis",
## reverting to default.
  • Fitting linear model: gamble ~ sex + status + income + verbal
      Estimate Std. Error t value Pr(>|t|)
    (Intercept) 7.631 12.93 0.5904 0.5582
    sexfemale -16.3 6.133 -2.657 0.01118
    status 0.1739 0.2083 0.8346 0.4088
    income 4.331 0.7636 5.672 1.265e-06
    verbal -1.802 1.614 -1.117 0.2707

La exclusión de la observación 24 genera varios cambios en todos los coeficientes del modelo. Ademas de ser un candidato a outlier es un valor muy influyente

  1. Check for influential points.
plot(model, 4)

halfnorm(cooks.distance(model), labs = teengamb_cl$.rownames)

El gráfico de distancias de Cook confirma que el valor 24 se aleja mucho del comportamiento de las otras variables, que es un valor influyente y probablemente un outlier.

influence(model)$coefficients %>%
    as_tibble(rownames = "id") %>%
    mutate(id = factor(id, levels = c(1:n))) %>%
    gather(key, value, -id) %>%
    ggplot(aes(fct_rev(id), value, group = key)) +
    geom_line(alpha = .75)+
    coord_flip() +
    geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
    facet_grid(~ key, scales = "free")

plot(model, 5)

pander(model)
Fitting linear model: gamble ~ sex + status + income + verbal
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 22.56 17.2 1.312 0.1968
sexfemale -22.12 8.211 -2.694 0.01011
status 0.05223 0.2811 0.1858 0.8535
income 4.962 1.025 4.839 1.792e-05
verbal -2.959 2.172 -1.362 0.1803 
teengamb_cl[24,] %>% pander()
  .rownames sex status income verbal gamble
24 24 male 27 10 4 156
  1. Check the structure of the relationship between the predictors and the response.

En primer lugar me remito al gráfico de dispersión de (a), que permite reflejar la distribución y la relación 1 a 1 de las variables.

car::avPlots(model)

En este caso el modelo refleja la relación directa del sexo masculino y del income con gamble, y la relación inversa suave con el score verbal. Es posible apreciar claramente que la observación 24 es un caso raro que se sale de las líneas de predicción para todos los coeficientes.

par(mfrow = c(2,2))
termplot(model, partial.resid = T)

En este caso los parciales de residuales reflejan lo que se ve en el primer gráfico de dispersión y es que hay variables que no tiene una distribución normal, como income y , que podrían hacer pensar en hacer que transformaciones de las variables mejorarían la estructura del modelo

3. (Ejercicio 3 cap. 6 pág. 97)

For the prostate data, fit a model with lpsa as the response and the other variables as predictors. Answer the questions posed in the first question.

n <- dim(prostate)[1]
p <- dim(prostate)[2]

prostate_cl <- 
    prostate %>%
    mutate(svi = as.logical(svi))

model <- 
    lm(lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45, data = prostate_cl)

ggpairs(prostate_cl, progress = F, lower = list(continuous = wrap("points", alpha = 0.5, size = .4)))

Lo primero que hago es crear una matriz de dispersión para tener una idea de como se distribuyen las variables y la relación 2x2 que tienen entre si. Utilizo la librería GGally que tiene la función ggpairs.

  1. Check the constant variance assumption for the errors.
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos

Gráficos Diagnósticos

Para evaluar la homocedasticidad del error utilizo los gráficos de:

  • residuales v/s valor ajustado
  • residuales estandarizados v/s valor ajustado.

En este caso no se observa un patrón definido por lo que se asume homocedasticidad.

  1. Check the normality assumption.
shapiro.test(resid(model)) %>% pander()
Shapiro-Wilk normality test: resid(model)
Test statistic P value
0.9911 0.7721

La gráfica Q-Q parece bastante pareja y aproximada a la linea. El test de Shapiro Wilk confirma la normalidad del error

  1. Check for large leverage points.
halfnorm(hatvalues(model), labs = rownames(prostate))

augment(model) %>%
    mutate(.rownames = rownames(prostate)) %>%
    filter(.hat > 2*(p/n)) %>%
    pander()
Table continues below
lpsa lcavol lweight age lbph svi lcp gleason pgg45 .fitted .se.fit .resid .hat
2.008 0.1823 6.108 65 1.705 FALSE -1.386 6 0 2.875 0.4072 -0.867 0.3305
2.158 1.423 3.657 73 -0.5798 FALSE 1.658 8 15 1.925 0.3311 0.2323 0.2184
2.298 0.6206 3.142 60 -1.386 FALSE -1.386 9 80 2.049 0.3478 0.2481 0.241
3.075 1.839 3.237 60 0.4383 TRUE 1.179 9 90 3.544 0.3098 -0.4689 0.1912
4.13 2.533 3.678 61 1.348 TRUE -1.386 7 15 4.07 0.3241 0.0593 0.2092
.sigma .cooksd .std.resid .rownames
0.7034 0.1227 -1.496 32
0.7119 0.004271 0.3709 37
0.7118 0.005703 0.402 41
0.7103 0.01423 -0.736 74
0.7124 0.0002605 0.09413 92

Para evaluar los valores extremos en el modelo, utilzamos dos gráficos.

  • La gráfica de los valores ajustados versus los residuales estandarizados (\(\hat{y} \times \sqrt{|{\hat{\epsilon}}|}\))
  • Una gráfica de media normal (halfnormal) versus los residuales estandarizados. De esta manera se pueden identificar los valores mas extremos en \(X\).

En este caso podemos ver que los valores para el caso 32, 37 y el 41 cumplen con las características para valores extremos.

  1. Check for outliers.
augment(model) %>%
    mutate(.t.resid = rstudent(model)) %>%
    filter(abs(.t.resid) > abs(qt(.05/(n * 2), n - p))) %>%
    pander()
Table continues below
lpsa lcavol lweight age lbph svi lcp gleason pgg45 .fitted
.se.fit .resid .hat .sigma .cooksd .std.resid .t.resid

Para valorar posibles outliers podemos utilizar la estrategia planteada en el libro de Faraway, y calcular los residuales studentizados. En este sentido los residuales cuyo valor superen un limite determinado por la corrección de Bonferroni, se podrían considerar posibles outliers.

En este caso no hay valores que cumplan con lo propuesto

  1. Check for influential points.
plot(model, 4)

halfnorm(cooks.distance(model), labs = rownames(prostate))

model_1 <- 
    prostate_cl %>%
    rownames_to_column() %>%
    filter(rowname != "32") %>%
    lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + lcp + 
    gleason + pgg45, data = .)

model_2 <- 
    prostate_cl %>%
    rownames_to_column() %>%
    filter(rowname != "47") %>%
    lm(formula = lpsa ~ lcavol + lweight + age + lbph + svi + lcp + 
    gleason + pgg45, data = .)
    


pander(model)
Fitting linear model: lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.6693 1.296 0.5163 0.6069
lcavol 0.587 0.08792 6.677 2.111e-09
lweight 0.4545 0.17 2.673 0.008955
age -0.01964 0.01117 -1.758 0.08229
lbph 0.1071 0.05845 1.832 0.0704
sviTRUE 0.7662 0.2443 3.136 0.002329
lcp -0.1055 0.09101 -1.159 0.2496
gleason 0.04514 0.1575 0.2867 0.775
pgg45 0.004525 0.004421 1.024 0.3089 
pander(model_1)
Fitting linear model: lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.1719 1.329 0.1293 0.8974
lcavol 0.5653 0.08847 6.39 7.93e-09
lweight 0.6217 0.202 3.077 0.002793
age -0.02127 0.01115 -1.908 0.05963
lbph 0.09559 0.05853 1.633 0.106
sviTRUE 0.7604 0.2426 3.135 0.002347
lcp -0.106 0.09036 -1.173 0.244
gleason 0.05069 0.1564 0.3241 0.7466
pgg45 0.004468 0.00439 1.018 0.3116 
pander(model_2)
Fitting linear model: lpsa ~ lcavol + lweight + age + lbph + svi + lcp + gleason + pgg45
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.0488 1.29 0.03783 0.9699
lcavol 0.5659 0.08615 6.568 3.574e-09
lweight 0.4577 0.1657 2.762 0.007009
age -0.01685 0.01095 -1.538 0.1277
lbph 0.1167 0.05711 2.043 0.04409
sviTRUE 0.8265 0.2394 3.452 0.0008625
lcp -0.09798 0.08876 -1.104 0.2727
gleason 0.1142 0.1562 0.7313 0.4665
pgg45 0.004457 0.004309 1.034 0.3038

Los casos mas influyentes son el 32 y el 47. Si sacamos del modelo la observación 32, el coeficiente de lweight cambia en \(.15\). Al quitar la observación 47, el modelo se preserva mejor.

influence(model)$coefficients %>%
    as_tibble(rownames = "id") %>%
    mutate(id = factor(id, levels = c(1:n))) %>%
    gather(key, value, -id) %>%
    ggplot(aes(fct_rev(id), value, group = key)) +
    geom_line(alpha = .75)+
    coord_flip() +
    geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
    facet_grid(~ key, scales = "free")

En este gráfico se ve como 32 genera un cambio importante en lweight. Llama la atención ademas la observación 69 que genera cambios en varios coeficientes y que tiene una distancia de cook importante.

plot(model, 5)

prostate_cl[c(69,32, 47),] %>% pander()
  lcavol lweight age lbph svi lcp gleason pgg45 lpsa
69 -0.4463 4.409 69 -1.386 FALSE -1.386 6 0 2.963
32 0.1823 6.108 65 1.705 FALSE -1.386 6 0 2.008
47 2.728 3.995 79 1.879 TRUE 2.657 9 100 2.569

Los casos 69, 32 y 47 son casos influyente. El caso 32 ademas tiene mucha palanca. Pero no tengo mucha evidencia como para decir que son outliers.

  1. Check the structure of the relationship between the predictors and the response.
car::avPlots(model)

La regresión parcial confirma que 32 y mas aún el69 son casos raros que se comportan de manera diferente a la mayoría en cuanto a las lineas de regresión parcial.

par(mfrow = c(2,2))
termplot(model, partial.resid = T)

Aquí destaca el caso 32 con un lweight muy alejado del resto, y la estrucutra de lbph (y en menor grado, de lcp, pgg45 y gleason) que no es para nada homogenea, lo que ya se veía en el grafico de dispersion.

4. (Ejercicio 4 cap. 6 pág. 97)

For the swiss data, fit a model with Fertility as the response and the other variables as predictors. Answer the questions posed in the first question.

n <- dim(swiss)[1]
p <- dim(swiss)[2]

model <- 
    lm(Fertility ~ Agriculture + Examination + Education + Catholic + Infant.Mortality, data = swiss)


swiss %>%
    ggpairs(progress = F, lower = list(continuous = wrap("points", alpha = 0.75, size = .8)))

Lo primero que hago es crear una matriz de dispersión para tener una idea de como se distribuyen las variables y la relación 2x2 que tienen entre si. Utilizo la librería GGally que tiene la función ggpairs.

  1. Check the constant variance assumption for the errors.
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos

Gráficos Diagnósticos

No se observa un patrón definido, por lo que se acepta una homocedasticidad del error.

  1. Check the normality assumption.
shapiro.test(resid(model)) %>% pander()
Shapiro-Wilk normality test: resid(model)
Test statistic P value
0.9889 0.9318

El gráfico *Q-Q se adapta bien a la linea. La prueba no permite rechazar la H0 de normalidad del error.

  1. Check for large leverage points.
halfnorm(hatvalues(model), labs = rownames(sat))

augment(model) %>%
    filter(.hat > 2*(p/n)) %>%
    pander()
Table continues below
.rownames Fertility Agriculture Examination Education Catholic Infant.Mortality .fitted .se.fit
La Vallee 54.3 15.2 31 20 2.15 10.8 50.74 4.246
V. De Geneve 35 1.2 37 53 42.34 18 34.8 4.838
.resid .hat .sigma .cooksd .std.resid
3.562 0.3512 7.221 0.03437 0.6172
0.2024 0.4558 7.254 0.0002047 0.03829

Para evaluar los valores extremos en el modelo, utilizamos dos gráficos.

  • La gráfica de los valores ajustados versus los residuales estandarizados (\(\hat{y} \times \sqrt{|{\hat{\epsilon}}|}\))
  • Una gráfica de media normal (halfnormal) versus los residuales estandarizados. De esta manera se pueden identificar los valores mas extremos en \(X\).

En este caso podemos ver que los valores para Vermont y Maine cumplen con las características para valores extremos. Llaman la atención también los puntos de La Valle y V. De Geneve que presentan valores extremos en una de las variables.

  1. Check for outliers.
augment(model) %>%
    mutate(.t.resid = rstudent(model)) %>%
    filter(abs(.t.resid) > abs(qt(.05/(n * 2), n - p))) %>%
    pander()
Table continues below
.rownames Fertility Agriculture Examination Education Catholic
Table continues below
Infant.Mortality .fitted .se.fit .resid .hat .sigma .cooksd .std.resid
.t.resid

No hay outliers por este criterio. Tampoco se observan en el gráfico de las distancias de Cook en los gráficos diagnósticos

  1. Check for influential points.
par(mfrow = c(2,2))
plot(model, 4)
halfnorm(cooks.distance(model), labs = rownames(swiss))

influence(model)$coefficients %>%
    as_tibble(rownames = "id") %>%
    gather(key, value, -id) %>%
    ggplot(aes(fct_rev(id), value, group = key)) +
    geom_line(alpha = .75)+
    coord_flip() +
    geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
    facet_grid(~ key, scales = "free") +
    scale_colour_viridis_d(end=.85) 

colMeans(swiss) %>% pander()
Fertility Agriculture Examination Education Catholic Infant.Mortality
70.14 50.66 16.49 10.98 41.14 19.94 
swiss %>%
    rownames_to_column() %>%
    filter(rowname %in% c("Porrentruy", "Sierre")) %>%
    pander()
rowname Fertility Agriculture Examination Education Catholic Infant.Mortality
Porrentruy 76.1 35.3 9 7 90.57 26.6
Sierre 92.2 84.6 3 3 99.46 16.3

Hay varias observaciones que resultan influyentes. POr le criterio de distancias de Cook, destaca Porrentruy que en segundo gráfico podemos ver, modifica agriculture y Infant.Mortality. Sierre modifica Infant.Mortality hacia el otro sentido, tiene una mortalidad baja, y Examination baja.

Hay dos puntos con alta palanca pero baja influencia, La Valle y V. De Geneve.

  1. Check the structure of the relationship between the predictors and the response.

En primer lugar me remito al gráfico de dispersión de (a), que permite reflejar la distribución y la relación 1 a 1 de las variables.

car::avPlots(model)

La estructura parece correcta. Podemos ver como Sierre se aleja de la tendencia de los modelos.

par(mfrow = c(2,2))
termplot(model, partial.resid = T, se = T)

Se pueden observar los valores extremos pero con poca influencia en Education y en Infant.Mortality Destaca la estructura bimodal de Catholic. Se podría realizar una transformación de la variable a una categorial o incluso binaria.

5. (Ejercicio 5 cap. 6 pág. 97)

Using the cheddar data, fit a model with taste as the response and the other three variables as predictors. Answer the questions posed in the first question.

n <- dim(cheddar)[1]
p <- dim(cheddar)[2]

model <- 
    lm(taste ~ Acetic + H2S + Lactic, data = cheddar)

cheddar %>%
    ggpairs(progress = F)

  1. Check the constant variance assumption for the errors.
par(mfrow = c(2,2))
plot(model)
Gráficos Diagnósticos

Gráficos Diagnósticos

No se ve un patrón el el gráfico de valores ajustados v/s residuales que haga pensar en problemas de varianza del error.

  1. Check the normality assumption.
shapiro.test(resid(model)) %>% pander()
Shapiro-Wilk normality test: resid(model)
Test statistic P value
0.9802 0.8312

El gráfico Q-Q sigue una distribución muy cercana a la recta y la prueba de S-W no permite desechar la H0 por lo que se puede asumir la normalidad de los residuales.

  1. Check for large leverage points.
halfnorm(hatvalues(model), labs = rownames(cheddar))

augment(model) %>%
    rowid_to_column() %>%
    filter(.hat > 2*(p/n)) %>%
    pander()
Table continues below
rowid taste Acetic H2S Lactic .fitted .se.fit .resid .hat
.sigma .cooksd .std.resid

Los casos 20 y 26 son los casos mas extremos.

  1. Check for outliers.
augment(model) %>%
    mutate(.t.resid = rstudent(model)) %>%
    filter(abs(.t.resid) > abs(qt(.05/(n * 2), n - p))) %>%
    pander()
Table continues below
taste Acetic H2S Lactic .fitted .se.fit .resid .hat .sigma
.cooksd .std.resid .t.resid

No hay outliers según este criterio.

  1. Check for influential points.
plot(model, 4)

halfnorm(cooks.distance(model), labs = rownames(cheddar))

influence(model)$coefficients %>%
    as_tibble(rownames = "id") %>%
    mutate(id = factor(id, levels = c(1:n))) %>%
    gather(key, value, -id) %>%
    ggplot(aes(fct_rev(id), value, group = key)) +
    geom_line(alpha = .75)+
    coord_flip() +
    geom_hline(yintercept = 0, linetype = 3, color = "darkred") +
    facet_grid(~ key, scales = "free")

Los casos 12 y 15 son casos influyentes.

model_1 <- 
    cheddar %>%
    rownames_to_column() %>%
    filter(rowname != "15") %>%
    lm(formula = taste ~ Acetic + H2S + Lactic, data = .)



pander(model)
Fitting linear model: taste ~ Acetic + H2S + Lactic
  Estimate Std. Error t value Pr(>|t|)
(Intercept) -28.88 19.74 -1.463 0.1554
Acetic 0.3277 4.46 0.07349 0.942
H2S 3.912 1.248 3.133 0.004247
Lactic 19.67 8.629 2.28 0.03108 
pander(model_1)
Fitting linear model: taste ~ Acetic + H2S + Lactic
  Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.76 17.62 -1.007 0.3233
Acetic -2.47 4.004 -0.617 0.5428
H2S 4.039 1.091 3.702 0.001061
Lactic 21.46 7.559 2.839 0.008862

El caso 15 modifica el coeficiente de Acetic en 2 , y Lactic en casi 2.

  1. Check the structure of the relationship between the predictors and the response.
car::avPlots(model)

par(mfrow = c(2,2))
termplot(model, partial.resid = T)

NO se ven mayores problemas estructurales en los gráficos.

6. (∗) (Ejercicio 6 cap. 6 pág. 98)

Using thehappydata, fit a model withhappyas the response and the other four variables aspredictors. Answer the questions posed in the first question.

7. (∗) (Ejercicio 7 cap. 6 pág. 98)

Using thetvdoctordata, fit a model withlifeas the response and the other two variables aspredictors. Answer the questions posed in the first question.

8. (∗) (Ejercicio 8 cap. 6 pág. 98)

For thedivusadata, fit a model withdivorceas the response and the other variables, except yearas predictors. Check for serial correlation.

9. (Ejercicio 3 cap. 7 pág. 110)

Using the divusa data:

  1. Fit a regression model with divorce as the response and unemployed, femlab, marriage, birth and military as predictors. Compute the condition numbers and interpret their meanings.
# Miramos colinearidad
divusa %>%
select(-year, -divorce) %>%
  cor() %>%
  corrplot.mixed(., lower = "number", upper = "circle", tl.pos = "lt")

model <- lm(divorce ~ unemployed + femlab + marriage + birth + military, data = divusa)

Podemos observar que si hay varios predictores que presenta correlación entre si.

X <- model.matrix(model)[,-1]
eig <- eigen(t(X) %*% X)

eig$values
## [1] 1174600.548   21261.741   16133.842    6206.181    1856.894
(c_nums <- sqrt(eig$values[1] / eig$values))
## [1]  1.000000  7.432684  8.532498 13.757290 25.150782

Los números de condición son relativamente pequeños (< 30) por lo que a pesar de la correlación de parejas encontrada, no parece existir un problema grave de colinearidad.

  1. For the same model, compute the VIFs. Is there evidence that collinearity causes some pre-dictors not to be significant? Explain.
VIF <- vif(X)
pander(rbind(vifs, SE = sqrt(vifs)))

Quitting from lines 789-791 (keymer_alejandro.Rmd) Error in rbind(vifs, SE = sqrt(vifs)) : object ‘vifs’ not found Calls: … withCallingHandlers -> withVisible -> eval -> eval -> pander -> rbind

Si bien hay algunas \(X_i\) donde la VIF se aleja de 1 que es la ortogonalidad, tampoco se aleja demasiado. Donde más se aleja es en \(X_{femlab}\) donde se puede interpretar que el error estándar aumenta en 1.9008618

model_1 <-
  lm(formula = divorce+2*rnorm(dim(divusa)[1]) ~ unemployed + femlab + marriage + birth + 
    military, data = divusa)


pander(model)
Fitting linear model: divorce ~ unemployed + femlab + marriage + birth + military
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.488 3.394 0.7331 0.4659
unemployed -0.1113 0.05592 -1.989 0.05052
femlab 0.3836 0.03059 12.54 1.106e-19
marriage 0.1187 0.02441 4.861 6.772e-06
birth -0.13 0.0156 -8.333 4.027e-12
military -0.02673 0.01425 -1.876 0.06471 
pander(model_1)
Fitting linear model: divorce + 2 * rnorm(dim(divusa)[1]) ~ unemployed + femlab + marriage + birth + military
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.956 5.229 0.374 0.7095
unemployed -0.06709 0.08616 -0.7786 0.4388
femlab 0.3795 0.04713 8.054 1.33e-11
marriage 0.1229 0.03762 3.267 0.001674
birth -0.1327 0.02403 -5.525 5.118e-07
military -6.257e-06 0.02195 -0.000285 0.9998

A pesar de eso, si agregamos un poco de ruido el modelo se mantiene relativamente estable.

  1. Does the removal of insignificant predictors from the model reduce the collinearity? Investigate.
model_1 <-
  lm(formula = divorce ~ unemployed + marriage + birth + military, data = divusa)

model_2 <-
  lm(formula = divorce ~ unemployed + femlab  + military, data = divusa)

# Modelo original
pander(model)
Fitting linear model: divorce ~ unemployed + femlab + marriage + birth + military
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 2.488 3.394 0.7331 0.4659
unemployed -0.1113 0.05592 -1.989 0.05052
femlab 0.3836 0.03059 12.54 1.106e-19
marriage 0.1187 0.02441 4.861 6.772e-06
birth -0.13 0.0156 -8.333 4.027e-12
military -0.02673 0.01425 -1.876 0.06471 
# Modelo 1, divorce ~ unemployed + marriage + birth + military
pander(model_1)
Fitting linear model: divorce ~ unemployed + marriage + birth + military
  Estimate Std. Error t value Pr(>|t|)
(Intercept) 42.01 2.246 18.7 2.369e-29
unemployed -0.5686 0.07551 -7.53 1.151e-10
marriage -0.05652 0.03566 -1.585 0.1173
birth -0.2289 0.02396 -9.552 1.978e-14
military -0.0155 0.02532 -0.6121 0.5424 
# VIF model 1
pander(vif(model_1))
unemployed marriage birth military
1.295 1.927 1.925 1.245 
# Modelo 2, divorce ~ unemployed + femlab  + military
pander(model_2)
Fitting linear model: divorce ~ unemployed + femlab + military
  Estimate Std. Error t value Pr(>|t|)
(Intercept) -3.784 1.232 -3.072 0.002988
unemployed 0.01835 0.06086 0.3015 0.7639
femlab 0.4412 0.02412 18.3 5.529e-29
military -0.008315 0.02014 -0.4128 0.681 
# VIF model 2
pander(vif(model_2))
unemployed femlab military
1.276 1.075 1.195

En este caso reviso dos modelos alternativos. Uno sin la variable femlab cuyo VIF era alta en l modelo original, y un segundo modelo con la variable femlab pero sin marriage ni birth que tienen una correlación de pareja alta con femlab y entre si.

Los dos modelos alternativos pierden bastante poder explicativo lo que se ve reflejado en un \(R^2\) mas pequeño, sobretodo el sin femlab. Por otra parte, es verdad que en el primer modelo sin femlab no baja la VIF (aumenta en las variables que tiene correlación, birth y marriage). En el segundo modelo, sin birth y marriage, la \(R^2\) baja, pero la VIF también.

Como conclusión, creo que el modelo original es correcto.

10. (Ejercicio 4 cap. 7 pág. 110)

For the longley data, fit a model with Employed as the response and the other variables as predictors.

  1. Compute and comment on the condition numbers.
model <- lm(Employed ~ ., data = longley)

X <- model.matrix(model)[,-1]
eig <- eigen(t(X) %*% X)

eig$values
## [1] 6.665299e+07 2.090730e+05 1.053550e+05 1.803976e+04 2.455730e+01
## [6] 2.015117e+00
(c_nums <- sqrt(eig$values[1] / eig$values))
## [1]    1.00000   17.85504   25.15256   60.78472 1647.47771 5751.21560

Los numeros de condición son bastante altos al igual que los eigenvalues, lo que sugiere multicolinearidad

  1. Compute and comment on the correlations between the predictors.
corrplot.mixed(cor(X), lower = "number", upper = "circle", tl.pos = "lt")

En este caso las variables predictoras tienen mucha mayor correlación entre si.

  1. Compute the variance inflation factors.
VIF <- vif(X)
pander(rbind(vifs, SE = sqrt(vifs)))

Quitting from lines 872-874 (keymer_alejandro.Rmd) Error in rbind(vifs, SE = sqrt(vifs)) : object ‘vifs’ not found Calls: … withCallingHandlers -> withVisible -> eval -> eval -> pander -> rbind

En este caso excepto por Armed.Forces, y algo menos Unemployed, la mayoría de las variables predictoras correlacionan entre si, y hacen que aumente el error estándar del modelo. Es probable que un modelo mas simple pueda funcionar mejor.

pander(model)
Fitting linear model: Employed ~ .
  Estimate Std. Error t value Pr(>|t|)
(Intercept) -3482 890.4 -3.911 0.00356
GNP.deflator 0.01506 0.08491 0.1774 0.8631
GNP -0.03582 0.03349 -1.07 0.3127
Unemployed -0.0202 0.004884 -4.136 0.002535
Armed.Forces -0.01033 0.002143 -4.822 0.0009444
Population -0.0511 0.2261 -0.2261 0.8262
Year 1.829 0.4555 4.016 0.003037

11. (Ejercicio 5 cap. 7 pág. 110)

For the prostate data, fit a model with lpsa as the response and the other variables as predictors.

model <- 
  lm(lpsa ~ ., data = prostate)
  1. Compute and comment on the condition numbers.
X <- model.matrix(model)[,-1]
eig <- eigen(t(X) %*% X)

eig$values
## [1] 4.790826e+05 6.190704e+04 2.109042e+02 1.756329e+02 6.479853e+01
## [6] 4.452379e+01 2.023914e+01 8.093145e+00
(c_nums <- sqrt(eig$values[1] / eig$values))
## [1]   1.00000   2.78186  47.66094  52.22787  85.98499 103.73114 153.85414
## [8] 243.30248

HAlgunos valores son altos y otros no, lo que sugiere que podría haber multicolinerarida en varias combinaciones lineales.

  1. Compute and comment on the correlations between the predictors.
corrplot.mixed(cor(X), lower = "number", upper = "circle", tl.pos = "d")

Se corrobora lo anterior en la medida de que hay alguno de los predictores que correlacionan bastante, como el caso de pgg45 con gleason y lcp; lcp con svl y lcavol.

  1. Compute the variance inflation factors.
VIF <- vif(X)
pander(rbind(vifs, SE = sqrt(vifs)))

Quitting from lines 919-921 (keymer_alejandro.Rmd) Error in rbind(vifs, SE = sqrt(vifs)) : object ‘vifs’ not found Calls: … withCallingHandlers -> withVisible -> eval -> eval -> pander -> rbind

De todas formas, a pesar de la correlación, no parece haber un mayor problema al analizar la VIF. La mayoría de los predictores suman poco error

12. (∗) (Ejercicio 8 cap. 7 pág. 111)

Use the fat data, fitting the model described in Section 4.2.

  1. Compute the condition numbers and variance inflation factors. Comment on the degree ofcollinearity observed in the data.
  2. Cases 39 and 42 are unusual. Refit the model without these two cases and recompute thecollinearity diagnostics. Comment on the differences observed from the full data fit.
  3. Fit a model withbrozekas the response and justage,weightandheightas predictors.Compute the collinearity diagnostics and compare to the full data fit.
  4. Compute a 95% prediction interval forbrozekfor the median values ofage,weightandheight.
  5. Compute a 95% prediction interval forbrozekforage=40,weight=200andheight=73. Howdoes the interval compare to the previous prediction?
  6. Compute a 95% prediction interval forbrozekforage=40,weight=130andheight=73. Arethe values of predictors unusual? Comment on how the interval compares to the previous twoanswers.

Ejercicios del libro de Carmona1.

(∗) (Ejercicio 9.1 del Capítulo 9 página 172)

Realizar el análisis completo de los residuos del modelo de regresión parabólico propuesto en lasección 1.2 con los datos de tráfico.2.

(∗) (Ejercicio 9.2 del Capítulo 9 página 172)

Realizar el análisis completo de los residuos de los modelos de regresión simple y parabólico pro-puestos en la sección 1.2 con los datos de tráfico, pero tomando como variable respuesta la velocidad(sin raíz cuadrada). Este análisis debe justificar la utilización de la raíz cuadrada de la velocidadcomo variable dependiente.